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HIGHLIGHTS 


►  We  derive  a  reduced-order  discrete-time  state-space  model  of  a  lithium-ion  cell. 

►  We  start  with  known  Laplace-domain  transfer  functions  of  the  porous-electrode  model. 

►  Then,  we  derive  transfer  functions  for  electrolyte  concentration  and  the  potentials. 

►  Using  the  discrete-time  realization  algorithm,  we  produce  a  reduced-order  model. 

►  A  fifth-order  model  can  accurately  model  a  dynamic  drive  cycle. 
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We  present  a  method  to  produce  a  physics-based  one-dimensional  discrete-time  state-space  reduced- 
order  model  (ROM)  of  a  lithium-ion  cell.  The  resulting  ROM  can  predict  the  five  variables  of  a  stan¬ 
dard  porous-electrode  model— reaction  flux,  solid  and  electrolyte  lithium  concentration,  and  solid  and 
electrolyte  potentials— at  any  location  across  the  cell  cross  section,  as  well  as  cell  terminal  voltage.  The 
method  to  generate  the  model  involves  first  linearizing  the  porous-electrode-model  equations,  and  then 
deriving  closed-form  Laplace-domain  transfer  functions  from  the  linearized  equations.  Next,  the 
discrete-time  realization  algorithm  (DRA)  is  used  to  convert  the  transfer  functions  into  an  optimal 
discrete-time  state-space  realization.  Advantages  of  this  approach  include  that  the  DRA  avoids  nonlinear 
optimization  and  gives  a  straightforward  method  for  selecting  the  system  order  for  the  ROM.  Simulation 
results  demonstrate  that  the  ROM  cell  voltage  predictions  and  the  ROM  internal  electrochemical  variable 
predictions  match  very  closely  with  results  obtained  by  simulating  the  full  nonlinear  porous-electrode 
partial  differential  equations. 

©  2012  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

An  important  requirement  for  a  battery  management  system  is 
that  it  be  able  to  compute  accurate  estimates  of  a  number  of 
fundamental  quantities  that  describe  the  present  status  of  the 
battery  pack’s  cells.  These  include:  state-of-charge  (SOC),  state-of- 
health  (usually  comprising  a  cell  resistance  and/or  capacity  esti¬ 
mate),  state-of-life,  available  energy,  available  power,  and  time-to- 
empty.  The  best  available  methods  to  compute  these  estimates 
require  high-fidelity  but  computationally  simple  mathematical 
models  of  battery  cell  input/output  (current/voltage)  dynamics. 
Furthermore,  we  believe  that  future  battery-management-system 
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requirements  will  necessitate  models  that  give  insight  into  cell 
internal  electrochemical  dynamics;  for  example,  to  predict  and  to 
minimize  degradation  mechanisms,  to  slow  aging. 

Broadly  speaking,  battery  models  that  might  be  used  in  such 
model-based  estimators  fall  into  two  categories:  (1)  empirical 
equivalent-circuit  type  models,  and  (2)  physics-based  models. 
Equivalent-circuit  models  (ECMs)  are  formulated  by  recognizing 
that  cell  current/voltage  behavior  is  often  well  approximated  by 
a  model  that  uses  a  voltage  source,  resistors,  capacitors,  and  War¬ 
burg  impedances  as  analogs  to  observed  behaviors  in  the  cell  [1], 
The  values  of  the  model  components  are  found  via  empirical 
system-identification  approaches,  based  on  lab-generated  cell-test 
data.  Present  battery-management-system  algorithms  rely  almost 
exclusively  on  ECMs  due  to  the  simplicity  and  robustness  of  these 
models,  which  allow  them  to  operate  in  real-time  embedded 
applications.  However,  ECMs  do  not  provide  insight  into  the  elec¬ 
trochemical  dynamics  that  occur  internal  to  the  cell  during 


J.L.  Lee  et  al.  /  Journal  of  Power  Sources  220  (2012)  430-448 


431 


operation,  and  this  insight  will  be  key  to  fully  utilizing  the  cell 
without  causing  premature  degradation. 

Alternately,  physics-based  models,  originally  developed  by 
Doyle,  Fuller,  and  Newman  [2,3],  give  equations  that  provide  full 
information  on  the  internal  electrochemical  dynamics  of  the  cell. 
However,  the  coupled  partial-differential  equations  (PDEs)  that 
form  the  model  result  in  computational  complexity  that  makes 
their  use  in  real-time  embedded  battery-management  applications 
impractical.  For  real-time  control  applications,  we  desire  the  best 
features  of  the  ECMs  and  the  physics-based  models;  namely,  we 
need  a  computationally  simple  and  robust  reduced-order  model 
(ROM)  of  cell  dynamics  that  predicts  not  only  the  cell  voltage  and 
SOC  of  the  cell  during  operation,  but  also  accurate  estimates  of  the 
cell’s  internal  electrochemical  variables. 

A  number  of  approaches  to  reduce  the  computational 
complexity  of  the  physics-based  models  have  been  explored  in 
the  literature.  The  primary  gains  have  been  made  by  observing 
that  much  of  the  computational  complexity  involved  in  solving 
the  rigorous  physics-based  model  equations  comes  from  solving 
for  the  concentration  of  lithium  in  the  solid  particles  of  the 
electrodes.  Many  methods  therefore  focus  on  making  approxi¬ 
mations  to  simplify  this  computation.  Both  Wang  et  al.  [4]  and 
Subramanian  et  al.  [5]  have  used  a  parabolic  profile  to  estimate 
the  concentration  throughout  the  solid  particle,  eliminating  the 
need  for  this  PDE.  This  approach  works  at  low  discharge  rates  but 
does  not  perform  well  for  highly  dynamic  current  inputs  such  as 
those  encountered  in  hybrid-  and  electric-vehicle  applications. 
Others  have  used  the  Pade  approximation  to  match  the  power 
series  representation  of  the  solid  diffusion  equation  to  a  desired 
order  [6],  The  Pade  approximation  matches  the  dc  values  but  also 
does  not  perform  as  well  at  higher  frequencies.  Cai  and  White 
use  proper  orthogonal  decomposition  [7]  to  find  the  electro¬ 
chemical  variables  at  discrete  locations  within  the  solid  particle 
and  across  the  ID  cell.  This  method  accurately  models  the  cell 
behavior  at  high  currents  but  requires  either  existing  experi¬ 
mental  data  or  simulation  results  to  generate  the  ROM.  As  cell 
parameters  change  due  to  aging,  we  desire  a  method  that  does 
not  require  offline  simulation  results  in  order  to  generate  an 
updated  ROM. 

In  our  opinion,  the  most  promising  approach  to  reduced-order 
modeling  has  been  introduced  by  Smith  et  al.  [8,9],  In  this  work, 
they  showed  how  to  linearize  the  coupled  PDEs  of  the  rigorous 
pseudo-2D  porous-electrode  model  and  derive  analytic  Laplace- 
domain  transfer  functions  from  the  linearized  model.  They  were 
able  to  model  several  internal  cell  electrochemical  variables  at 
different  spatial  locations  within  the  cell  this  way:  reaction  flux 
j(z,t),  solid-electrolyte  potential  difference  </>s_e(z,t),  overpotential 
77 (z,t),  and  solid  particle  surface  concentration  cs,e(z,t).  They  did  not 
derive  an  analytic  transfer  function  for  the  concentration  of  lithium 
in  the  electrolyte  ce(x,t);  as  an  alternative,  they  used  the  finite 
element  method  to  solve  for  this  variable. 

Smith  et  al.  then  used  a  nonlinear  optimization  routine  to  fit 
a  reduced-order  model  to  these  transfer  functions  in  such  a  way  as 
to  minimize  the  error  between  the  frequency-domain  response  of 
the  reduced-order  model  and  the  full-order  linear  transfer-function 
model.  In  principle,  this  nonlinear  optimization  step  can  work  well. 
In  practice,  however,  it  often  does  not.  The  results  of  nonlinear 
optimization  are  very  sensitive  to  the  initial  starting  estimate,  they 
are  not  guaranteed  to  be  globally  optimal,  are  not  guaranteed  to 
converge  in  bounded  time,  and  the  model  order  must  be  found  via 
trial  and  error.  Nonlinear  optimization  is  not  well  suited  for  unsu¬ 
pervised  operation  in  embedded  systems,  which  limits  its  appli¬ 
cability  in  fielded  battery-management  applications  where  the 
parameters  of  the  transfer  functions,  and  hence  the  ROM  itself,  may 
need  to  adapt  as  the  cell  ages. 


In  this  paper,  we  propose  an  improvement  to  this  approach.  We 
begin,  as  did  Smith  and  colleagues,  by  finding  transfer  functions  of 
the  cell  variables  of  interest.  However,  we  remove  the  problematic 
nonlinear  optimization  step  and  instead  find  the  reduced-order 
model  using  the  discrete-time  realization  algorithm  (DRA)  [10], 
which  overcomes  the  major  disadvantages  listed  above.  Namely, 
the  DRA  computes  a  discrete-time  state-space  model  directly  from 
transfer  functions,  without  requiring  an  initial  guess  of  the 
reduced-order  model  parameters,  while  yielding  a  globally  optimal 
result  in  bounded  time,  and  also  supplying  guidance  to  choosing 
the  order  of  the  reduced-order  model.  In  [10],  Lee  et  al.  demon¬ 
strate  the  DRA  on  a  single-input  single-output  system.  Here,  we 
show  how  to  use  the  method  to  form  a  single-input  multiple- 
output  cell  model.  Additionally,  we  extend  the  model  of  Smith  et  al. 
by  deriving  analytic  transfer  functions  for  the  solid  potential,  0s(z,t), 
electrolyte  potential,  0e(x,t),  and  ce(x,t).  We  also  modify  their 
derivation  steps  somewhat  to  allow  for  a  particle  surface-film 
resistance  Rfiim  in  the  Butler-Volmer  equation,  and  we  propose 
nonlinear  output  corrections  that  improve  the  accuracy  over 
a  strictly  linear  model. 

This  paper  is  the  second  in  a  series  of  planned  works.  The  first 
paper  in  the  series  introduced  the  DRA  as  a  method  for  converting 
a  transcendental  transfer  function  into  an  optimal  discrete-time 
state-space  reduced-order  model  [10],  This  is  the  second  paper  in 
the  series,  which  shows  how  to  find  transcendental  transfer  func¬ 
tions  corresponding  to  lithium-ion  internal  cell  dynamics.  The 
scope  of  the  work  in  this  paper  is  limited  to  showing  how  to  create 
a  ROM  that  is  linearized  to  give  accurate  predictions  around  a  single 
SOC  and  temperature  set-point  only.  A  planned  third  paper  in  this 
series  will  show  how  to  enhance  the  method  presented  herein  to 
allow  wide  variations  in  temperature  and  SOC. 

This  paper  is  organized  in  the  following  way.  In  Section  2,  we 
reproduce  the  five  equations  of  the  standard  pseudo-2D  porous- 
electrode  model.  We  then  derive  transfer  functions  for  the  elec¬ 
trode  variables  j(z,t),  r)(z,t),  cSie(z,t),  and  <j>s(z,t )  using  a  methodology 
similar  to  that  in  [8,9],  but  which  is  modified  to  allow  the  addition 
of  a  particle  surface-film  resistance.  We  then  show  how  to  model 
the  electrolyte  variables  (j>e(x,t)  and  ce(x,t)  as  transfer  functions.  In 
Section  3,  we  briefly  review  the  DRA  and  show  how  to  use  it  to  yield 
optimal  reduced-order  discrete-time  state-space  realization  of 
these  transfer  functions.  In  Section  4,  we  describe  how  the  model’s 
voltage  estimates  can  be  improved  over  purely  linear  estimates  by 
adding  nonlinear  output  corrections.  Results  are  presented  in 
Section  5  where  we  compare  predictions  from  our  reduced-order 
model  to  those  computed  using  the  PDEs  of  the  rigorous  pseudo- 
2D  model. 

2.  Derivation  of  transfer  functions  of  reduced-order  model 

In  this  section,  we  derive  transfer  functions  corresponding  to 
linearized  equations  of  the  internal  electrochemical  dynamics  of 
a  cell.  The  derivation  is  quite  long;  however,  the  final  reduced-order 
model  produced  using  these  transfer  functions  is  computationally 
very  simple — similar  to  the  complexity  of  an  ECM.  Furthermore,  the 
transfer  functions  are  capable  of  very  accurate  predictions  of 
internal  cell  electrochemical  states  and  also  cell  terminal  voltage. 

2.1.  Pseudo-2D  porous-electrode  model  of  cell 

We  begin  by  presenting  the  partial  differential  equations, 
boundary  conditions  and  initial  conditions  that  comprise  the 
pseudo-2D  porous-electrode  cell  model  developed  by  Doyle,  Fuller 
and  Newman  [2,3],  which  is  the  starting  point  to  our  derivation.  We 
reproduce  these  equations  because  everything  that  follows 
depends  on  their  exact  form,  and  variations  in  how  these  equations 
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are  presented  in  the  literature  might  otherwise  introduce  an 
ambiguity  that  results  in  substantial  confusion. 

The  pseudo-2D  porous-electrode  cell  model  is  written  in  terms 
of  the  electrochemical  variables  that  summarize  the  internal  state 
of  a  cell  during  operation.  Specifically,  they  model  electrode-scale 
variables:  reaction  flux  j,  potential  in  the  solid  phase  </>s,  and 
concentration  of  lithium  in  the  solid  phase  cs.  They  also  model  cell- 
scale  variables:  potential  in  the  electrolyte  phase  <pe  and  concen¬ 
tration  of  lithium  in  the  electrolyte  phase  ce. 

Solid  particles  in  the  electrodes  are  modeled  as  spheres  with 
radius  Rs.  The  concentration  of  lithium  in  these  particles  cs(r,x,t)  is 
assumed  to  be  radially  symmetric  within  any  given  particle,  so 
overall  is  a  function  of  the  radial  distance  r  from  the  center  of  the 
particle  as  well  as  the  spatial  location  x  of  the  particle  along  the 
one-dimensional  cross-section  across  the  cell  (x  =  0  at  the  negative 
current  collector  and  x  =  L  at  the  positive  current  collector).  Lithium 
in  the  solid  is  assumed  to  move  via  diffusion: 

dcs(r,x,t)  _  Ds  6  fadcs(r,x,t)\ 

at  —  r2  ar  V  dr  )'  {  > 


The  boundary  conditions  to  this  PDE  are 
<UQ.0>  and  d*Cs(Rs,x4 


where  the  reaction  flux  j(x,t)  is  a  measure  of  the  amount  of  lithium 
moving  across  the  boundary  of  the  solid  particle.  The  initial 
concentration  profile  of  lithium  in  the  particle  is 

cs(r,x,0)  =  cSi0,  0  <r<Rs. 

Concentration  of  lithium  in  the  electrolyte  ce(x,t)  is  derived  from 
the  principle  of  mass  conservation,  and  is  found  to  be 


where  the  boundary  conditions  are  given  by 

ace(Q,t)  9Ce(L,t) 

ax  ax 


The  initial  concentration  of  lithium  in  the  electrolyte  is  ce(x,0)  =  ce,o 
for  0  <  x<L. 

Conservation  of  charge  in  the  solid  phase  gives  the  equation  for 
the  potential  in  the  solid  phase  0s(x,t),  which  is 

£  (<76ff^s(X’ t})  “  aMX’ t]  =  °’  (3) 
with  boundary  conditions 
„eff80s(O ,t)  _  I 

ax "  '  ^  a  ax  A  ~  Iapp' 

Conservation  of  charge  in  the  electrolyte  phase  gives  the  equa¬ 
tion  for  the  potential  of  the  electrolyte  phase  </>e(x,t),  which  is 

£  (Keff^0e(*’t})  4  (*0^  "*(*.0)  +asFj(x,t)  =  0,  (4) 

where  Keff  =  xcemg  and 


4ff  = 


dln/±\ 
dlncey ’ 


with  the  boundary  conditions 

30e(O,t)  a  0e(l,t) 

ax  ax 

In  this  work  we  assume  the  value  for  dln/±/dlnce  is  set  to  0. 

Finally,  the  Butler— Volmer  equation  gives  the  reaction  flux  j(x,t) 
at  different  points  in  the  electrodes 

j(X,t)  =  kcl  “(cmax-Cs,e)'~ac“e  X  (exp((1 

-exp(‘S”))'  <5) 

where  the  overpotential  t)  is 

n  =  <t>s-<t>e-  Uocp  -JfRflim.  (6) 

Eqs.  (1)— (5)  comprise  the  pseudo-2D  porous-electrode  model. 
The  main  dimension  is  the  cross-sectional  distance  x  across  the  cell; 
the  pseudo  dimension  is  the  radial  distance  r  into  a  particular 
electrode  particle.  Similar  to  Smith  et  al.,  we  eliminate  the  need  for 
computing  the  pseudo  dimension  r  by  using  transfer-function 
methods.  We  then  show  how  to  compute  transfer  functions  for 
the  cell  variables  at  any  desired  set  of  x  locations  across  the  cell. 


2.2.  Transfer  functions  for  reaction  fluxj(z,t),  overpotential  rj(z,t), 
and  solid  surface  concentration  cs  e(z,t) 

We  begin  by  finding  transfer  functions  for  three  cell  variables 
that  are  very  tightly  coupled  together:  the  reaction  flux  j(z,t),  the 
overpotential  r](z,t),  and  the  solid  surface  concentration  cs,e(z,t). 
Note  the  introduction  of  the  dimensionless  spatial  variable  z,  which 
is  defined  only  within  the  electrodes  and  not  within  the  separator. 
In  the  negative  electrode,  z  =  x/Ln  and  in  the  positive  electrode, 
z  =  (I-x)/Ip.  In  both  cases,  z  —  0  at  the  outer  current  collector/ 
electrode  interface  and  z  =  1  at  the  inner  electrode/separator 
interface.  This  is  illustrated  in  Fig.  1,  which  shows  the  set  of  variables 
that  exist  in  each  domain  of  the  cell’s  cross  section,  and  whether  the 
variable’s  position  variable  is  generally  chosen  to  be  x  or  z. 

We  embark  on  the  derivation  of  these  three  linear  transfer 
functions  from  the  porous-electrode  equations  in  a  manner  largely 
similar  to  that  of  Smith  et  al.  [8,9],  although  we  approach  the 
linearization  differently  in  order  to  add  a  particle  surface  film 
resistance  term  f?fiim  to  the  Butler— Volmer  equation.  We  make  the 
same  fundamental  two  assumptions  as  [8,9]:  (1)  we  assume  linear 
behavior  and  enforce  that  linear  behavior  by  linearizing  nonlinear 
expressions  as  needed,  and  (2)  we  assume  that  reaction  flux  can  be 
decoupled  from  electrolyte  concentration,  which  is  equivalent  to 
assuming  that  electrolyte  potential  is  independent  of  electrolyte 
concentration  when  deriving  the  transfer  function  for  reaction  flux. 


Negative  electrode  ( n )  Separator  (m)  Positive  electrode  (p) 


0j(z,  0  1) 

Mx,t) 

Ux,  t)  Uz,  t) 

Cs.e{z,  t)  Ce(x,  t) 

Ce(x,t ) 

ce(x,t)  cs.e(z,i) 

jn(z,  t) 

) 

n  fa  0 

nfa  0 

0 

t) 

z=0  z=  1 

z  =  1  z  =  0 

x  =  0  x=Ln  x  =  Ln+Lm  x=L 

Fig.  1.  Schematic  diagram  of  cell  showing  variables  that  exist  in  each  domain  and  the 
orientation  of  z  in  the  electrodes. 
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We  start  with  Assumption  1,  and  linearize  the  Butler— Volmer 
equation  around  the  set-point  p*  defined  as 

P*  =  {cs,e  =  CSi0,Ce  =  Ce> o,0s_e  =  Uocp  (cS,o)J  =  0}, 


cs,e(z,  t)  =  CS(RS,Z;  t)  =  Cs(Rs,Z,  t)  -  cs,0.  (12) 

Eq.  (1),  modified  to  be  written  in  terms  of  the  debiased 
variable,  is 


where  we  have  defined  0s_e  =  0s-0e-  Re-arranging  Eq.  (5)  gives 


Approximating  the  left-hand-side  (LHS)  of  Eq.  (7)  with  the  first 
two  terms  of  its  Taylor-series  expansion  gives 


9/ 

=  o  +  o  X  (cs,e-cs,o)  +  0  x  (ce  —  ce  0) 


=  j/jo, 


kclfia (frnax  -  Cs  0) 1  “c“0 


where  jo  =  o“(cmax  -  cs  0)'_ac"0.  Similarly,  we  expand  77  in  the 
right-hand-side  (RHS)  of  Eq.  (7)  using  Eq.  (6),  and  approximate  it 
with  the  first  two  terms  of  its  Taylor-series  expansion  to  get 


RHS  =  exp((1  RTa)f  (0s_e  -  Uocp  -  AW)) 

-eXP(-RT  ~~  Uocp  “ 

=  RHs(p*)  +  ^  (cs,e  -  cSi0) 

9RHSI  r.  ,  ,,  9RHSI  . 

+  ^S_e "  Uocp(Cs.°))  +~dT\p.J 

=  0  +  RT  ^s_e  “  Uocp(Cs.°)) 


_  F  |9t/0c 
RT  hlCse 


(Cs,e  -  Cs,o) 


We  now  remove  the  constant  bias  from  0s_e  and  cs  e  by  defining 
debiased  versions  as  0s_e  =  0s_e  -  Uocp(cSjo)  and  cs,e  =  cSie  -  cSjo, 
where  the  tilde  decoration  is  used  throughout  this  paper  to 
refer  to  a  debiased  variable.  We  substitute  these  debiased  variables 
into  our  results,  and  set  Eq.  (8)  equal  to  Eq.  (9)  to  arrive  at  the 
linearized  Butler— Volmer  equation 


l 

jo 


F  rat/ocpi  1 
Rr[icSie  y 


Cs,e  - 


F2Rfj|m. 
RT  J' 


Rearranging  Eq.  (10)  to  solve  for  tj>s_e  gives 


(10) 


UDll(r2 

r2  9r  V 


with  boundary  and  initial  conditions 


-j(z,t), 


Cs(r,z,0)  =  0,  0  <  r  <  Rs. 

A  transfer-function  solution  to  this  partial  differential  equation 
was  derived  by  Jacobsen  and  West  [11],  and  is 


CSle(s)  Rsf  tanh(/3)  \ 
J(s)  Ds  \tanh((3)  -  fi)  ’ 


(13) 


where  /3  =  Rs  ^s/Ds. 

We  will  now  use  the  two  results  found  to  this  point,  Eqs.  (11) 
and  (13),  to  derive  a  transfer  function  for  0s_e.  We  initially  focus 
on  the  negative  electrode  only,  and  later  generalize  to  the  positive 
electrode  as  well.  In  the  normalized  z  dimension,  the  solid  potential 
equation  (Eq.  (3))  is  written  as 

=  a*Fj'  (M) 


with  boundary  conditions  ((<7efr/Ln)(90s/5z))|  q  =  -iapp/A  and 
(90s/0z)  [z=1  =  0.  The  potential  in  the  electrolyte  is 

Keff  92 

=  -QsFj>  (15) 

where  Assumption  2  is  used  to  omit  the  electrolyte-concentration 
dependent  term  from  the  left-hand  side  of  Eq.  (4)  and  where 
a  linearized  /ceff  is  defined  as  K(ce,o)£emE  >n  the  transfer  function 
models  only.  The  boundary  conditions  are  (90e(z,  t)/9z)|z=0  =  0 
and  (Keff/In)(90e(z!  t) /3z)  |  —  iapp/A.  Subtracting  Eq.  (15)  from 

Eq.  (14)  gives  a  single  equation  for  the  phase  potential  difference 

0s-e 


|U-.  =  (’6) 

with  boundary  conditions  ((<Teff/In)(00s_e/0z))|  =  ((-Keff/Ln) 

(00s_e/0z))  =  -iapp/A.  We  take  the  Laplace  transform  of 

Eq.  (11)  7=1 


<Ps-e  =  F(Rct  +  RfilmU  + 


$S-e(Z,S)  =  FRtot](Z,S)  +  J 


CSle(Z,S), 


(17) 


where  we  define  the  charge-transfer  resistance  as  Rct  —  RTI(joF2). 
We  can  further  condense  notation  by  defining  Rtot— Rct  +  Rfiim  to  get 

0s_e(Z,  t)  =  FRtotj(Z,  t)  +  J  Cs,e(Z,  t).  (11) 


We  leave  this  result  to  the  side  for  the  time  being  and  now  focus 
on  the  solid  surface  concentration.  We  define  cSie(z,t)  =  cs(Rs,z,t), 
and  corresponding  debiased  versions 


and  then  rewrite  the  rightmost  term  as 


From  Eq.  (13),  we  have  a  transfer  function  for  CSje(s)/J(s)  at 
a  single  spatial  location  in  the  electrode.  However,  note  that 
diffusion  happens  only  along  the  radial  dimension  r  and  not 
the  lateral^  dimension  z,  so  we  can  write  the  same  transfer 
function  Cs,e(s)/J(s)  at  every  z  location.  In  other  words, 
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Cs,e{z,s)/J(z,s)  =  Cs,e{s)/J(s)  for  every  value  of  z.  Inserting  this 
result  into  Eq.  (17),  we  get 


J%p> 


S  (  tanh(ff) 

I  FDS  \tanh(/?)-/3 


Now,  writing  Eq.  (16)  in  the  Laplace  domain  gives 


with  boundary  conditions, 

(7eff  d&s-e(Z,  s)  I  -Keff  90s  e(2,  S)  I  -Japp (s) 

Ln  0Z  Ln  9Z  A 

lz=0  lz=l 


For  convenience  of  notation,  we  define  dimensionless  variable 
*s)  as1 


Ks) 


Ln\/^ 


Qs 

Keff 


[atJocpI  1  (  tanh(ff)  \ 

[  ^  IcJ  FD^  Vtanh(^)  B) 


(19) 


which  allows  us  to  write 

a20s_e(z,s)  9,  ~  .  , 

-V(?)4>s-e(Z,S)  =  0. 

The  solution  to  this  homogeneous  differential  equation,  after 
satisfying  both  boundary  conditions,  is 

0s_e(Z,S)  Ln  /  1  ,  ,  ,  ,  , 

Japp (s)  “  Al/(s)sinh(Ks))  l,«effC0Sh(i'(S)Z) 

+  Jfrcosh(Ks)(z-l))j.  (20) 

From  this  result,  we  find  the  transfer  function  for  the  reaction 
flux  to  be 

J(Z,S)  =  J(2,S)  0S-e(Z,^'.  #|s)  #s-e(z,s) 

,apP(S)  0s_e(z,  s)  ,apP(s)  asFL2(^-L  +  -Lj  ,app(s)  ’ 

where  J(z,s)/0s-e(z,s)  is  found  from  Eq.  (18).  Expanding  gives 

Jfes)  K*) 

Japp  (s)  asFLnA(Keff  + 

( <refrcosh(i/(s)z)  +  Keffcosh(i'(s)(z  -  1))\  , 

\  sinh  (v(s))  )■  121 J 


This  is  the  first  transfer  function  we  have  set  out  to  find,  and  is 
equivalent  to  the  corresponding  result  in  [8,9]  except  for  the  Rflim 
term  embedded  inside  of  our  p(s). 

The  overpotential  rj  is  given  by  Eq.  (6)  and  repeated  below, 

V  =  0s  —  0e  -  JJocp  —  jFRfiim  ■ 

We  use  Eq.  (11)  and  the  definition  of  0s_e  —  0S  -  <Pe  — 

0 s-e  +  Uocp(CSio)  to  get 

0s-e  =  FRtotj  +  JJocp (cs,o)  +  |  QC°CP |  Jcs.e, 
and 

7 1  =  FRcd  ~  JJocp  (Cs,e)  +  JJocp (cs,o)  +  |  qc  P|  J  Cs,e- 


We  note  that  J/0cp(cs,e)  ~JJocp(cSi0)  +  [(aJ/0cp/9cs,e)|^ps,e,  leaving 
us  with  the  following  transfer  function  for  rj. 


V{z,s) 
Japp  (s) 


=  FRct 


J(z,s) 

Japp(s) 


(22) 


Finally,  the  transfer  function  for  the  surface  concentration  is  given 
by 

Cs,e(Z,S).  Cs,e(Z,S)  J(Z,S) 

Japp(S)  J(Z,  S)  Japp(s) 

v(s)Rs  (  tanh(/3)  \ 

asFLnADs(n^  +  efm)  \tanh(/3)  -  Bj 
<7effcosh(p(s)z)  +  Keffcosh(y(s)(z  -  1)) 
sinh(p(s)) 

Eqs.  (21)— (23)  are  the  first  three  transfer  functions  used  in  the 
reduced-order  cell  model.  These  transfer  functions  are  valid  for  the 
negative  electrode.  For  the  positive  electrode,  the  derivation 
follows  exactly  the  same  steps,  and  the  resulting  transfer  functions 
are  the  same  as  those  for  the  negative  electrode  region,  but  are 
multiplied  by  -1,  and  the  negative-electrode  parameters  are 
replaced  by  their  corresponding  positive-electrode  values. 


2.3.  Transfer  function  for  solid  potential  <t>s(z,t) 


We  now  proceed  to  derive  an  independent  transfer  function  for  the 
solid  potential  0s(z,t)  in  the  negative  electrode.  From  Eq.  (14)  we  have 


o^_&_ 

J2  8z2 


0s 


=  asFj. 


Integrating  both  sides  of  the  equation  from  0  to  z  and  using  the 
boundary  condition -((<reff/Ln)(a/az)0s(z,t))|  q  =  iapp/A  gives 


o 


Taking  the  Laplace  transform  and  dividing  both  sides  of  the  equa¬ 
tion  by  Japp(s)  gives, 


1  We  keep  the  notation  of  Smith  et  al.  and  use  the  Greek  letter  “nu”  for  this 
variable,  but  always  write  it  as  i/(s)  to  distinguish  it  from  cell  voltage,  which  is 
written  as  v(t). 


a  Ss(z,s)  _  -Ln  L2asF  fM,s)  r 
azJapp(s)  Acreff  +  o^ff  J  Japp(s)  v 
0 


We  define  a  debiased  variable  0s(z,  t)  =  <j>s(z,  t)  -0S(O,  t),  and 
integrate  both  sides  of  Eq.  (24)  to  give 
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$s (z,s)  _  Ln/ceff(cosh(j'(s))  -cosh((z-  l)v(s))) 
4pp (s)  ~  Affefr(f£eff  +•  ffefr)p{s)sinh(j'(s)) 

Ln<reff(l  -  cosh(z?'(s))  +23'(s)sinh(j'(s))) 
+  ff=ff)j'(s)sinh(p(s)) 


(25) 


In  order  to  solve  the  PDE  for  <j>e  we  integrate  Eq.  (26)  to  arrive  at 

*e(x,t)-*e(0,t)  = /  (Z^F+2P(1-f°)  9lnCaf’°)  dg-  (27) 
o 


As  before,  the  positive-electrode  version  of  this  transfer  function  is 
simply  multiplied  by  —  1,  and  negative-electrode  parameters  are 
replaced  by  their  corresponding  positive-electrode  values. 


2.4.  Transfer  function  for  electrolyte  potential  cj)e(x,t) 


We  now  develop  an  independent  transfer  function  for  0e  as 
a  function  of  x,  the  non-normalized  length.  Since  we  are  dealing 
with  variables  that  exist  over  the  entire  cell  width  and  not  over 
a  single  electrode  only,  we  also  introduce  subscripts  “n”  to  denote 
parameter  values  corresponding  to  the  negative  electrode,  “m”  to 
denote  values  in  the  separator  region,  and  “p”  to  denote  values  in 
the  positive  electrode,  whenever  making  a  distinction  is  necessary.2 
Starting  with  Eq.  (4),  and  integrating  once  with  respect  to  x  gives 


9</>e(x,t)  _  —ie(x,  t)  2RT  /  0\  0lnce(x,t) 

0X  T*  ,+  fV  L+)  0X  ’ 


(26) 


where 


J  asFjtf,  t)d£j,  0  <  x  <  In; 

4pp(t)/A  Ln  <  X  <  Ln  +  Lm; 

'.mV-  j  asFM,t) df,  Ln  +  Lm<x<L. 


In  the  rigorous  electrode  model,  the  effective  ionic  conductivity,  xeff,  is 
a  function  of  ce.  For  the  reduced-order  model,  /cefr  is  constant  and 
determined  from  the  equilibrium  ce  value.  The  transfer  function  for  the 
ionic  current  is  found  by  integrating  the  previously  derived  equation 
for  the  reaction  flux  transfer  function.  In  the  negative  electrode  this  is 


Ie(X,S) 

4pp(s) 


asF 


J 

o 


4pp(s) 


df 


^sinh^-^sinh^-^^) 
A(xeff  +  <7eff)  sinh(j'(s)) 

+  A(xeff  +  <reff)' 


The  transfer  function  for  ie(x,t)  in  the  separator  is  /e(x,s)//app(s)  =  1/ 
A.  The  transfer  function  for  the  positive  electrode  is 


4(x,s) 

4pp(s) 


o*sinh(<L-*M)  +^sinh(^+L-- 
A(/ceff  +  ffeff)sinh(?'(s)) 

Keff 

+  A(xeff  +  freff)- 


xMsy 


2  The  subscript  “m"  is  used  to  denote  the  electrolyte  phase  in  the  separator  in  the 
"middle”  of  the  cell,  as  "s”  is  already  used  to  denote  the  solid  phase. 


We  define  <j)e(x,  t)  =  0e(x,  t )  -  0e( 0,  t).  0e(x,  t)  is  comprised  of  two 
parts  corresponding  to  the  two  parts  of  the  equation  above: 

fcM,  -  /  "i^°« 

o 

The  first  part  [0e(x,  t)]j  is  solved  using  transfer  functions.  The 
second  part,  [0e(x,  t)]2  will  be  determined  from  the  value  of  ce(x,t), 
which  we  find  in  Section  2.5.  The  transfer  function  for  the  first  part, 
for  locations  x  in  the  negative  electrode,  is 

[*»H,  -i  ’[km* 

fapp(s)  Knff  J  ^app(s) 

0 


Ln(^j  i1  ~  cosh(X^n(S)))  -  wn(s)sinh(j'n(s)> 
A^ff  +  o*ff)  pn(s)sinh(%(s)) 

4i  f  COSh(l'n(s) )  -  COShf^”  ~ ^Vn^\\ 

+  — - 7 - c - V - n - y/ .  (28) 

A(xf +  <rf)pn(s)sinh(iin(s)) 

For  locations  x  in  the  separator  region,  the  transfer  function  is 

(29) 

4PP(s)  A<f +  A(xeff+<)pn(s) 

Finally,  for  locations  x  in  the  positive  electrode,  the  transfer  func- 


Vanlt^11^)  -  fn(s)  j 


4pp(s) 


A(xeff  +  ffefr)yn(s) 

+^fcosh(*'p(s))) 

A(x|ff  +  <tpff )  sinh  (t'p  (s) )vp(s) 

ipCosh((k±^4W) 

A(xf  +  <7|ff)sinh(j'p(s))!'p(s) 

,  V<  ^  Lp  n  , 

A(xeff  +  sinh(pp(s))i'p(s)  A^f  +  <7eff) 
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For  the  second  term  of  Eq.  (27)  we  have 


The  homogeneous  problem  is  written  as 


o  (31) 

2RT(1  ~tP+)  (ce{x,t)\ 

F  m\ce(0,t)J- 

In  our  final  reduced-order  model,  we  will  be  interested  in  values  of 
the  electrolyte  potential  at  particular  locations  across  the  cell.  At 
each  x  location,  we  will  then  have  a  transfer  function  to  solve  for 
[®e(x,  s)]  T .  This  transfer  function  will  be  from  Eq.  (28),  Eq.  (29)  or  Eq. 
(30)  depending  on  whether  x  is  in  the  negative  electrode,  separator, 
or  positive  electrode,  respectively.  In  the  next  section,  we  derive  an 
analytic  transfer  function  for  the  electrolyte  concentration  ce(x,t), 
which  we  need  to  be  able  to  compute  [</>e(x ,  t)]2. 

2.5.  Transfer  function  for  the  electrolyte  concentration  ce(x,t) 


0CeM  =  (D|ff(x)  9Ceax C))  •  (33) 

In  addition  to  the  boundary  conditions  given  above  we  have 
internal  boundary  conditions  where  the  three  regions  of  the  cell 
join.  These  internal  boundary  conditions  are 

Ce((Ln  +  Lm)-,t)  =  Ce((Ln  +  Im)+,t) 

r  3c-M  r 

n  ax  ax 

D  pac.((L„  +  L„)M) 

m  ax  P  ’$fk. 


In  [8,9],  Smith  et  al.  used  transfer  functions  to  solve  for  j,  (J>s_e,  V, 
and  cs,e,  but  then  used  the  finite-element  method  to  solve  for  ce.  We 
prefer  to  use  the  transfer  function  approach  for  all  variables,  since 
this  lets  us  jointly  optimize  the  reduced-order  model  over  all 
transfer  functions  using  the  DRA.  So,  in  this  section,  we  focus  on 
deriving  a  transfer  function  for  ce(x,t). 

The  governing  equation  for  the  electrolyte  concentration  is 
given  by  Eq.  (2),  repeated  here  but  now  explicitly  noting  the 
x-dependence  of  re  and  D|ff, 

te(x) dCe^’  ^  (of  (X) ^  Ce(X,  t) j  +  (f  - 1° ) a<f(x:  t ),  (32) 

with  boundary  conditions 


Using  the  method  of  separation  of  variables,  the  solution  to  Eq. 
(33)  can  be  split  into  a  function  of  time  only  and  a  function  of 
position  only 

Ce(X,t)  =  h(t)V(x). 

Substituting  this  assumed  form  into  Eq.  (33)  gives 


dh(t) 

dt 


W(x)  = 


1  0 
te(x)  0X 


Df(x)h(t) 


d  W(x) 


Separating  the  time-dependent  variables  on  one  side  of  the  equa¬ 
tion,  the  position  dependent  variables  on  the  other,  we  get 


^  =  0  and  ^  =  0, 

0X  0X 

and  uniform  initial  condition  ce(x,0)  =  ce,o- 

The  analytic  transfer  function  for  the  electrolyte  concentration  is 
derived  using  the  following  steps.  First,  we  find  the  homogeneous 
solution  to  Eq.  (32)  using  the  method  of  separation  of  variables  and 
by  finding  eigenfunctions  that  are  orthonormal  with  respect  to  the 
weighting  function  re(x).  The  solution  to  the  inhomogeneous 
problem  is  found  by  projecting  the  concentration  function  onto 
these  eigenfunctions  to  find  the  generalized  Fourier  coefficients. 
Finally,  we  derive  the  Laplace  domain  transfer  function  that  will 
solve  for  the  Fourier  coefficients. 


1  dh(t) 
h(t)  dt 


1  d 

ce(x)W(x)  dx 


(of(x) 


d> m\ 

dx  )' 


Since  one  side  of  the  equation  is  a  function  of  time  only,  and  the 
other  side  is  a  function  of  position  only,  and  the  equations  are  equal 
to  each  other  for  all  values  of  t  and  x,  then  both  sides  must  be 
constant: 


1  dh(t) 
h(t )  dt 


1  d 

ce(x)V(x)  dx 


(nf(x) 


dW(x)\ 

ry 


=  -A. 


Splitting  this  into  separate  equations  in  time  and  position  gives 


2.5.1.  Solution  to  the  homogeneous  problem 

We  first  use  the  method  of  separation  of  variables  to  find  the 
solution  to  the  homogeneous  problem.  Our  approach  is  similar  to 
[12]  but  whereas  Schmidt  et  al.  required  D|ff  and  i:e  to  be  constant 
across  the  entire  cell,  we  allow  for  different  values  in  each  of  the 
three  regions  of  the  cell.  In  order  to  simplify  notation,  we  use 
(cn,  X  <  Ln; 

£e(x)  =  <  £m,  Ln<X<Ln  +  Im; 

[cp,  Ln  +  Lm<X<L. 

The  values  of  rn,  £m,  and  ep  are  constant.  Likewise,  the  electrolyte 
diffusivity  is  assumed  to  be  constant  within  each  region  and  we 
simplify  the  notation  using 

f  Dn,  X  <  Ln; 

Df(x)  =  Dm,  Ln  <  X  <  Ln  +  Lm; 

[Dp,  Ln  +  Lm  <  X  <  L. 


d-fl  -  -m  (34) 

$  (“'"W  Trj"  (35) 

There  are  infinitely  many  specific  eigenvalues  A  that  satisfy  these 
equations,  with  specific  corresponding  eigenfunctions  W.  We 
change  notation  so  that  h(t)  >->h(t;  A)  and  change  W(x)^>W(x-,  A).  The 
solution  to  Eq.  (34)  is  then  of  the  form 

h(t;A)  =  Zj(0;  A)e-At. 

Eigenfunction  solutions  to  Eq.  (35)  have  three  distinct  parts:  one 
each  for  the  negative  electrode,  the  separator,  and  the  positive 
electrode.  For  the  negative  electrode,  we  have 
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The  boundary  condition  at  x  =  0  eliminates  the  sine  term,  which 
leaves 


which  solve  the  homogeneous  PDE.  We  denote  specific  solutions  as 
7k,  and  search  for  solutions  as  roots  of 


-  =  -k5sin(  J*h)  +k6c  os(  J^t]  =  0.  (37) 


Wn(X;X)  =  kj  cos 

In  the  separator  and  positive  electrode,  the  eigenfunctions  are 

Wm(x-,X)  =  k3cosU^xj  ^x) 

Wp(x-,X)  =  ^)+Wn(g. 


Using  the  internal  boundary  conditions,  we  can  solve  for  k3  and 
k4  in  terms  of  k\  by  solving  the  linear  equations 


In  a  similar  manner,  we  solve  for  ks  and  fe6  in  terms  of  k3  and  k 4 
where  x*  =  Ln  +  Lm: 


It  is  not  possible  to  find  a  closed-form  solution  for  7k  as  ks  and  k6  are 
themselves  functions  of  7k.  Instead,  we  do  a  numerical  search  of  the 
zero  crossings  of  Eq.  (37)  until  we  find  for  the  desired  number  of 
eigenvalues.3  We  denote  the  ordered  set  of  eigenvalues  as  (7k). 
Then,  the  solution  to  the  homogeneous  problem  is  the 
superposition 

Ce(x,t)  =  h( 0;7k)¥(x;7k)e-Akt. 

k= 0 


2.5.2.  Solution  to  the  inhomogeneous  problem 

We  now  generalize  the  homogeneous  solution  to  the  inhomo¬ 
geneous  case  to  solve 


dce(x,  t) 
at 


l  a 

£e(x)  dX 


(of(x) 


a Ce(X,  t)\ 
8X  J 


(38) 


We  transform  ce(x,t)  into  a  series  expansion  using  the  eigen¬ 
functions  from  the  previous  section  as  a  basis  set 


Ce(X,t)  =  ^C*ik(t)lF(X;7k), 
fc=0 


(39) 


where  c*  k(t)  are  the  generalized  Fourier  coefficients  of  ce(x,t).  Taking 
the  partial  derivative  of  this  equation  with  respect  to  time  gives 


3 Ce(X,  t)  _  y'  dCe,k(t) 


<F(X;7k). 


Substituting  Eq.  (40)  into  Eq.  (38)  gives 


(40) 


The  eigenfunction  over  the  entire  cell  cell  is  then 

(*Pn(X;7),  0  <  X  <  Ln; 

*frm(X; 7),  Ln  <  X  <  Ln  +  Lm;  (36) 

IFp(X;7),  in  +  Lm  <  X  <  L. 

By  Sturm-Liouville  theory  [13],  the  different  eigenfunctions  are 
automatically  orthogonal  with  respect  to  the  weighting  function 
£e(x),  and  k]  is  chosen  to  satisfy 


£%tW;4) 


To  reduce  the  left-hand  side  of  the  above  equation,  we  can 
multiply  both  sides  by  W(x;km)ee(x)  and  integrate  from  0  to  L. 


r  “  dc* .  (t) 

J  E  £  V(X;  AkW{X]  X™^e(x)dx 


0 

L 

+  /as(l-  ttyj(x, t)W(x; 7m)dx. 


J  W2(x;  7)re(x)dx  =  t, 
o 

making  them  furthermore  orthonormal  with  respect  to  te(x).  The 
final  boundary  condition  31fr(L;7)/5x  =  0  allows  us  to  find  the  set  of  7 


3  Note  that  one  aspect  of  Sturm-Liouville  theory  guarantees  that  the  kth 
eigenfunction  >P'(x;Ak)  has  k  zero  crossings  in  (0,L),  so  by  plotting  the  eigenfunctions 
over  the  cell  width  it  is  simple  to  detect  if  the  numeric  search  has  succeeded  (cf. 
Fig.  4  for  an  example  of  what  the  eigenfunctions  can  look  like).  In  practice,  we  find 
that  relatively  few— between  about  five  and  ten— eigenvalue/eigenfunction  pairs 
are  needed  to  approximate  ce  very  well. 
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The  eigenfunctions  are  orthonormal  with  respect  to  £e(x)»  so  the 
integral  on  the  left-hand  side  is  non-zero  only  when  m  =  n.  This  gives 

o 

L 

+  /*(!_  t°)j(x,  t)W(x;  Ak)dx.  (41 ) 

o 

Using  Green’s  identity  [13]  on  the  first  term  on  the  right-hand 
side  of  the  equation  gives, 


where  we  have  defined 


4)  =  J  as(\-tl)j(x,t)W(x-,Xk)dx. 


2.5.3.  Transfer  function 

To  compute  ce(x,t)  via  transfer-function  methods,  we  note  that 
ce(x,t)  is  formed  from  a  summation  of  c*  k(t)  terms.  We  can  compute 
ce(x,t )  if  we  know  c*k(t),  so  we  proceed  by  finding  a  transfer 
function  for  c*  k(t): 


(42) 

The  right-hand  side  goes  to  zero  because  of  the  boundary  conditions 

Mir(x;Afc)|  =  8ce(*,t)|  =  Q_ 

dx  lxe{0,L}  dx  lxe{0,L} 

The  left-hand  side  of  Eq.  (42)  has  two  parts.  The  first  part  can  be 
written  as 

/*<*■'> 

°  (43) 

=  -Afc  J  ce(x,t)W(x;Xk)ee(x) dx, 

0 

because  If^xiAfc)  satisfies  the  homogeneous  case  Eq.  (35) 


dc4(t)  *  * 

dt  -4ce,k(t)  +h(t ) 

sCe,k(S)!  lS~4Ceik(s)  +Jk(s) 

^,k(s)  i  rk(s) 

fapp  (S)  S  +  Ak  fapp  (S) 

We  further  examine  jk(t), 

L 

fk(t)  =  J  as  (l  -  t+^)j(x,  t)W(x-,  Ak)dx 
0 
In 

*  J  Qs(l  -  t+)j(x,tW(X;^k)dX 
0 

L 

+  J  as  (l  -  t°)j(x,  t)W(x;  Ak)dx 

I„+t m 


+AkEe(x)W(X;Ak)  =  0. 

Since  the  right-hand  side  of  Eq.  (42)  is  zero,  we  can  write 

J  W(x;  Ak)  ^  ^Dgff(x)  °Ce^’ C)  j  dx 

0 

I 

=  -Ak  J  Ce(X,t)£e()i)y(X;  Ak)dx. 

0 

Using  this  in  Eq.  (41 )  gives 

dc?k(t)  r 

— —  =  -4  J  ce(x,  t)W(x:  Ak)re(x)dx 
o 

i 

+  J  as(l  -t°)j(x,t)»P(x;Ak)dx 

0 

L 

=  -4 4,k(t)  +  J  as(l-t°)j(x,t)1P(x;Ak)dx 


=  +Jk°S*(t)- 


For  the  negative  electrode, 


JkegV)  =  Os(l  -  t°)  Ji(x,t)¥(x;Ak) dx 

0 


(44) 


Note  that /(z,s)//app(s)  is  given  by  Eq.  (21 ),  where  z  =  x/Ln  in  the 
negative  electrode.  Performing  the  integration  of  Eq.  (44)  gives, 


Jkeg*(s)  4  (1  -  f?)InSin(Ln)  (<ff +  ^ffCOSh(j'n(s)))pn(s) 
Japp(S)  ^.eff  +  ffeff^  ^L*)2+p2(s)  j  Sinh(Pn(s)) 

ki  (1  -  t°)  ^K®ff-|-(7^ffCOS(L*)  jkn(s) 

AF(Kf  +  <f)((L*)2+p2(s)) 


~4ce,k(^)  +ik(t;  4)> 
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where  L*  =  Ln\AnAk/Dn.  For  the  positive  electrode, 


2.6.  Summary  of  transfer  Junctions  forming  the  model 


fkos'(t)  =  as(l  -t°)  J  j(x,t)T(x;4) dx 

tn+tm 


Jr'® 

/app(S) 


/ 


J{(L-x)/Lp,s) 

/app(S) 


W(x;  Ak)dx, 


Before  proceeding,  we  quickly  summarize  the  transfer  func¬ 
tions  of  the  model.  Recall  that  z  is  a  normalized  unitless  spatial 
variable  in  the  electrodes  that  takes  on  value  0  at  the  current 
collectors  and  1  at  the  electrode/separator  interface,  and  that  x  is 
the  spatial  variable  across  the  cell  and  takes  on  value  0  at  the 
negative  electrode  current  collector  and  L  at  the  positive  electrode 
current  collector.  Common  to  all  of  the  transfer  functions  is 
Eq.  (19),  which  gives  the  dimensionless  function  v(s)  in  terms  of 
0  =  Rss/Ws- 


which  gives, 


;P°s’(s)  _Ml  -  pCOs(I*)  (af  +  Kfcosh(vp(s))yp(s) 
faPP(s)  AF^tcf  +  ^L*)2+t/2(s)jsinh(rp(s)) 

M1  -  t°)LpSin(Lp)  (xf  +  <j^frcosh(pp(s)))i'p(s) 
AF(xf  +  of'j  (  (l;)  2+^(s))  sinh(j/p  (s)) 

+  ke(l  -t?)L  pCos(L*m)  (xf  +  ofcosh(vp(s))yp(s) 

AF(xf  +  (  (l;) 2+vj  (s))  sinh(pp(s)) 

ks  (1  -  t+)L*p  (sin  (I*)  (of  +  nfcosh(vp(s)))  )  vp(s) 
AF(xf  +  ((l;) 2+^(s))sinh(i/p(s)) 

k5(l  - t£)  (<j£ffcos(l*m)  +K|ffcos(L*))^(s) 

AF(<f  +  ffeff)((L;)2+p2(s)) 

k6(l  - 1+)  (ffpffsin(LJm)  +  Kpffsin(L*)^p(s) 

- - - 2 - ' - ,  (46) 

AF(K|ff  +  <rf)((L*)  +J-2(s)) 

where  L*  =  Lpf  cpkk/Dp,  L*m  =  (L„  +  Lm)  ^  rp),k/Dp,  and 
L*  =  l^J £pAk/Dp.  The  transfer  function  for  the  generalized  Fourier 
coefficients  of  the  electrolyte  concentration  is 


c;k(s)  j 

4pp(5) 

Since  the  concentration  is  given  by  Eq.  (39)  we  can  find 
a  transfer  function  for  ce  at  any  x  location  across  the  cell.  This 
transfer  function  is  given  by 


kes\s)  |  r>)l 

(4pp(s)  fapp(s)  J 


Ce(X,S) 
4pp  (S) 


to'aPPW 


(48) 


where  JC  is  the  number  of  eigenvalue/eigenfunction  expansion 
terms  chosen  to  approximate  the  exact  solution.  The  first  term  in 
the  summation,  (C*0(s)//app(S))|f,'(x;  Ak),  is  equal  to  the  steady-state 
concentration  of  lithium  in  the  electrolyte  phase,  ce,o.  We  define 
ce(x.  t)  =  ce(x.  t)  -  ce  o  and  remove  this  term  from  the  summation 
to  give 


Ce(X,S) 

lapp(s) 


f.c;k(s) 


*P(X;4). 


(49) 


The  following  transfer  functions  model  cell  parameters  at  any 
desired  sets  of  spatial  locations  in  the  cell  based  only  on  the  cell 
input  current  /app(s). 

Reaction  flux:  The  local  reaction  flux  transfer  function  J(z,s)/ 
/app(s)  in  the  negative  electrode  is  given  by  Eq.  (21 ).  In  the  positive 
electrode,  simply  multiply  this  transfer  function  by  -1,  and  replace 
negative-electrode  parameters  in  the  equation  by  the  corre¬ 
sponding  positive-electrode  quantities. 

Overpotential:  The  local  overpotential  transfer  function  ?;(z,s)/ 
fapp(s)  in  the  negative  electrode  is  given  by  Eq.  (22).  In  the  positive 
electrode,  simply  multiply  this  transfer  function  by  -1,  and  replace 
negative-electrode  parameters  in  the  equation  by  the  corre¬ 
sponding  positive-electrode  quantities. 

Solid  surface  concentration:  Lithium  concentration  in  the 
solid  is  denoted  as  cs(rg,t).  We  are  concerned  here  only  with  the 
concentration  of  lithium  at  the  surface  of  the  particles,  as  that  is 
what  determines  reaction  rate.  The  surface  concentration  is 
denoted  as  cSie(z,t).  Furthermore,  we  define  a  debiased  surface 
concentration  that  subtracts  out  the  equilibrium  concentration: 
cs,e(z,t)  =  cs,e(z,t)-cs,o.  The  transfer  function  of  the  debiased 
surface  concentration  Cs,e(z,s)//app(s)  in  the  negative  electrode  is 
given  by  Eq.  (23).  In  the  positive  electrode,  simply  multiply  this 
transfer  function  by  —  1,  and  replace  negative-electrode  parame¬ 
ters  in  the  equation  by  the  corresponding  positive-electrode 
quantities.  And,  once  the  debiased  surface  concentration  is 
computed,  the  true  surface  concentration  can  be  found  as 
c*,e(z,  t)  =  cs,e(z,  t )  +  cs-0. 

Potential  in  solid:  Solid  potential  in  an  electrode  is  denoted  as 
(j>s(z,t).  We  define  a  debiased  solid  potential  that  subtracts  out  the 
potential  at  the  electrode’s  current  collector  where  z  =  0: 
<j>s(z,t)  —  <t>s{z,t)  -  0S (0,  t).  The  transfer  function  $s(z,s)//app(s)  of 
the  debiased  solid  potential  in  the  negative  electrode  is  given  by  Eq. 
(25).  In  the  positive  electrode,  simply  multiply  this  transfer  func¬ 
tion  by  —  1,  and  replace  negative-electrode  parameters  in  the 
equation  by  the  corresponding  positive-electrode  quantities.  And, 
since  we  define  the  potential  of  the  negative-electrode  current 
collector  to  be  zero,  we  have  0s(z,t)  =  0s(z,t)  in  the  negative 
electrode.  In  the  positive  electrode,  we  have  that 
0s(z,  t)  =  0s(z,  t)  +v(t),  where  v(t)  is  the  cell  voltage  and  can  be 
found  as  described  in  Section  4.2. 

Potential  in  electrolyte:  Electrolyte  potential  is  denoted  as 
0e(x,t).  We  define  a  debiased  electrolyte  potential  that  subtracts  out 
the  potential  at  the  negative  current  collector  where  x  =  0: 
0e(x,  t)  —  0e(x,  t)  —  0e(O,  t).  Furthermore,  we  break  up  0e(x,  t)  into 
two  parts:  <j>e{x,t)  —  [0e(x.  t)],  +  [0e(x.  t)]2.  The  transfer  function 
[<£e(x,s)]i//app(s)  for  x  locations  in  the  negative  electrode  is  given  by 
Eq.  (28).  For  x  locations  in  the  separator,  this  transfer  function  is 
given  by  Eq.  (29),  and  for  x  locations  in  the  positive  electrode,  this 
transfer  function  is  given  by  Eq.  (30).  Once  the  appropriate 
[0e(x;  t)]  ]  is  computed,  we  compute  <Pe(x,  t)  =  [0e(x.  t)] , 
+(2RT(1  -  £®)/F)ln(Ce(x,t)/Ce(0,t)),  where  ce(x,t)  is  computed  as 
described  in  the  “concentration  in  electrolyte”  summary  item. 
Finally, 
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0e(X,  t)  =  <j)e(X,  t)  +  </>e(0,  t) 

=  0e(X,t)  -  ^s_e(0,  t), 


where  0s_e( 0,  t)  atx  =  0  is  found  via  the  transfer  function  in  Eq.  (20) 
using  negative-electrode  variables. 

Concentration  in  electrolyte:  Lithium  concentration  in  the 
electrolyte  is  denoted  ce(x,t).  We  use  an  eigenfunction  expansion  to 
express  ce(x,t)  in  terms  of  K+l  weighted  partial  transfer  functions, 
where  K  is  at  least  as  large  as  the  anticipated  reduced-order-model 
system  order.  The  first  term  of  the  expansion  is  the  equilibrium 
concentration  so  we  define  ce(x,  t )  =  ce(x.  t )  —  ce,o  and  use  the 
transfer  function  in  Eq.  (49)  to  find  Ce(x,s)//aPP(s).  Details  of 
finding  the  eigenfunctions  W(x;Ak)  are  given  in  Section  2.5.  The 
individual  partial  transfer  functions  are  expressed  in  Eq.  (47).  The 
two  parts  of  this  transfer  function  are  found  in  Eqs.  (45)  and  (46). 


Note  that  the  transfer  function  that  is  being  converted  from 
a  continuous-time  representation  to  a  reduced-order  discrete-time 
state-space  realization  may  be  single-input  single-output,  multi¬ 
input  single-output,  single-input  multi-output,  or  multi-input 
multi-output.  Here,  we  are  considering  the  case  of  a  single  input 
iapp(t)  and  multiple  outputs,  where  the  outputs  comprise  the  set  of 
transfer  functions  from  Section  2.2  that  the  user  would  like  to 
implement.  The  overall  transfer  function  being  implemented  by  the 
DRA  is  then  the  vertical  concatenation  of  all  sub-transfer  functions. 
For  example,  if  the  user  would  like  to  determine  j(0,t)  and  77(1,1),  the 
overall  transfer  function  would  be 


Ms) 

fapp(S) 

V(1,s) 
hpp  (s) 


3.  Discrete-time  realization  algorithm 

We  have  now  derived  all  of  the  necessary  continuous-time 
transfer  functions  required  to  model  the  internal  dynamics  of 
a  lithium-ion  cell.  Proceeding,  our  goal  is  to  convert  these  transfer 
functions  into  an  optimal  reduced-order  discrete-time  state-space 
model  of  the  form 

x(t  +  Ts)  =  Ax(t)  +  Biapp(t) 
y(t)  =  Cx(t)  +  Dlapp(t), 

where  Ts  is  the  sampling  period  of  the  discrete-time  ROM,  x(f)  is  the 
“state”  vector  of  the  model  at  time  t,  y(t)  is  the  “output"  vector  of 
the  model  at  time  t,  and  A,  B,  C,  and  D  are  matrices.  Note  that  the 
ROM  states  and  outputs  are  valid  only  at  values  of  time  that  are 
integer  multiples  of  Ts. 

The  challenge  in  accomplishing  this  conversion  is  that  the 
transfer  functions  are  transcendental  in  the  Laplace  variable  s, 
which  makes  closed-form  solutions  intractable.  To  overcome  this 
obstacle,  we  introduced  the  “discrete-time  realization  algorithm” 
(DRA)  in  the  first  paper  in  this  series  [10].  The  DRA  relies  on  the 
Ho— Kalman  algorithm  [14],  which  generates  an  optimal  reduced- 
order  discrete-time  state-space  realization  from  a  discrete-time 
system  pulse  response.  Therefore,  the  DRA  first  finds  the  discrete¬ 
time  pulse  response,  and  then  uses  the  Ho-Kalman  algorithm  to 
compute  the  state-space  realization.  A  summary  of  the  steps 
implemented  by  the  DRA  is  given  below.  A  detailed  description  of 
the  algorithm  is  given  in  [10], 

1  Sample  the  continuous-time  transfer  function  in  the  frequency 
domain  at  a  high  rate  Fs  (not  related  to  Ts),  and  take  the  inverse 
discrete  Fourier  transform  (1DFT)  to  get  an  approximation  to 
the  continuous-time  impulse  response. 

2  Form  the  continuous-time  step  response  from  the  continuous¬ 
time  impulse  response,  and  then  compute  the  discrete-time 
pulse  response  values  from  the  continuous-time  step 
response,  assuming  a  sample  and  hold  circuit  connected  to  the 
system  input. 

3  Generate  a  discrete-time  state-space  realization  using  the 
deterministic  Ho— Kalman  algorithm.  This  algorithm  returns 
the  reduced-order  A,  B,  and  C  matrices  from  the  discrete-time 
pulse-response  sequence  in  Step  2.  The  order  of  the  system  is 
determined  from  the  sorted  singular  values  of  the  Hankel 
matrix  that  is  computed  as  part  of  the  algorithm.  The  D  matrix 
is  found  by  the  initial  value  theorem. 

4  Transform  the  state  space  system  into  the  desired  final  form 
using  a  similarity  transformation,  if  required. 


Prior  to  implementing  DRA  Step  1,  H(s)  must  be  strictly  stable 
(e.g.,  it  may  not  contain  a  pole  at  s  =  0,  which  corresponds  to 
integration  dynamics).  If  H(s)  has  a  pole  at  s  —  0,  this  pole  must  first 
be  removed.  To  do  so,  we  determine  the  residue  associated  with 
this  pole  using  res0  =  lims  >0sH(s).  Then,  we  run  the  DRA  on  the 
strictly  stable  system  H*(s)  where  H*(s)  =  H(s)-reso/s.  For  the 
transfer  functions  in  Section  2,  two  have  a  pole  at  the  origin:  the 
potential  difference  transfer  function  (Eq.  (20))  and  the  solid 
surface  concentration  transfer  function  (Eq.  (23)).  We  calculate  the 
residue  for  these  poles  and  then  subtract  them  out  giving, 


— e(Z,S)  _  ^-e 
fapp(s)  4pp(S)  S 


(50) 


Cs,e(Z,S)  =  Cs,e(Z,S) 
4pp(s)  *app(s)  S 

where 

_3  [^ocpl 

2res0  [  9cs,e  cs>0 

<?s-e  asFIARs 


(51) 


(52) 


-3 

asFLARs ' 


(53) 


The  superscript  asterisk  on  the  transfer  function  denotes  that 
the  transfer  function  had  a  pole  at  the  origin  that  is  now  removed. 
We  use  the  transfer  functions  in  Eq.  (50)  and  Eq.  (51)  instead  of  Eq. 
(20)  and  Eq.  (23)  as  input  to  the  DRA.  In  the  final  state  space 
representation,  we  augment  the  state-space  realization  with  an 
integrator  state  to  account  for  removal  of  the  pole  in  the  DRA.  The 
reduced-order  approximation  to  the  original  system  via  the 
augmented  discrete-time  state-space  model  is: 


Adra 

0  1 

[  0 

1  J 

[cDra 

reso  J 

"  4PP(0 


(54) 


(55) 


where  reso  is  a  column  vector  containing  the  residues  of  the 
transfer  functions  with  a  pole  at  the  origin.  The  subscript  “DRA”  on 
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the  A,  B,  and  C  matrices  denotes  the  state  space  realization 
computed  by  Step  3  of  the  DRA. 

The  length  of  the  output  vector  depends  on  the  number  of  z 
locations  (for  the  electrode-scale  properties)  and  x  locations  (for 
the  cell-scale  properties)  chosen  by  the  user  for  evaluation.  For 
example,  solving  for  the  seven  variables  0s_e,  j,  77,  cse,  0S,  [<Pe]  1  and 
ce  at  four  spatial  locations  each  (e.g.,  both  current  collectors  and 
both  electrode-separator  interfaces)  produces  a  total  of  28  outputs. 
The  output  y(t)  then  has  the  following  structure: 


y  (t) 


4>S-e(Z,  t) 
j(Z,t) 
V(z,t) 
Cs,e(Z,t) 
fs(Z,t) 

[Mx,  t)]  ^ 

Ce(X,t) 


where  each  variable  is  a  four-vector  corresponding  to  the  four 
spatial  locations  evaluated. 

Note  that  DRA  Step  1  requires  we  select  a  high-rate  sampling 
frequency  Fs  to  approximate  the  impulse  response  of  the 
continuous-time  system.  We  find  that  the  results  are  not  very 
sensitive  to  Fs,  but  that  it  is  important  that  the  duration  of  the 
truncated  impulse  response  is  long  enough  to  capture  the  slow 
time  constants  in  the  model. 

While  the  elements  of  the  output  vector  y(t)  of  the  DRA  have 
physical  meaning,  the  only  state  variable  within  the  state  vector 
x(t)  that  has  independent  physical  interpretation  is  the  integrator 
state.  If  this  state  is  Xi(t),  then  we  have  that  the  average  concen¬ 
tration  of  lithium  in  the  solid  is  equal  to 


for  the  negative  electrode  (with  a  similar  result  for  the  positive 
electrode).  If  we  make  the  simplifying  assumption  that  reaction 
flux  is  approximately  uniform  across  the  electrode,  then  the 
“correct”  nonlinear  average  can  be  computed  using 
j(t)  —  i3pp(t) / (asFL„A)  and  by  inverting  the  Butler— Volmer  equa¬ 
tion  to  find  tjnoniinear-  Substituting  j(t)  for  j(t)  in  Eq.  (5)  gives 


f_aF-  \\ 

P I  gj  ^nonlinear  I  J  > 


where  j0(z,  t)  =  kcl~a(cm3x  -  cs)l_“cf  with  values  of  ce  and  cs 
corresponding  to  this  particular  cell  location  at  this  point  in  time.  If 
a  =  0.5,  the  “correct”  nonlinear  average  overpotential  at  a  point  in 
the  negative  electrode  is 


T/nonlinear(Z'  0  — 


2RT  (  iapp(t)  \ 
IF  \2j0(z,t)asFLnAj- 


The  nonlinear  corrected  rj  is  computed  by  subtracting  7?iinear(f) 
from  the  linear  DRA  output  ))DRA(z,t),  and  then  by  adding 

’Inonlinear^,  t) 

V{Z,t)  =  JJdRa(zi  0  —  ^linearM  +  ^nonlinear!2)  £)•  (56) 


Solid  surface  concentration:  The  debiased  solid  surface 
concentration  variable  cSie(z,  t)  is  defined  as  cs,e(z,  t)  = 
cs,e(z,  t)  -  cSio  (Eq.  (12)).  Therefore,  the  value  for  cs,e(z,t)  is  found  by 
adding  the  equilibrium  concentration  to  cs,e(z,  t): 


Cs,e(Z,t)  =  Cs,e(Z,t)  +  Csfi. 


Cs,avg(t)  -  Cs, 0  +  (c£f  ^(t). 

As  the  average  concentration  of  lithium  in  an  electrode  is  related  to 
the  state-of-charge  of  the  electrode  via  an  affine  transformation,  the 
value  of  the  integrator  state  is  key  to  being  able  to  estimate  cell  SOC. 

4.  Nonlinear  output  equations 

The  reduced-order  model  produced  by  the  DRA  is  entirely  linear, 
whereas  true  cell  behavior  is  nonlinear.  In  some  cases,  the 
nonlinear  cell  variable  is  simply  an  affine  function  of  the  linear 
output,  and  the  nonlinear  value  is  recovered  by  adding  a  constant  to 
the  linear  output.  In  other  cases,  the  transformation  is  more 
complex.  This  section  first  shows  how  to  recover  the  nonlinear 
internal  cell  variables  from  the  linear  cell  variables,  then  how  to 
estimate  the  nonlinear  cell  voltage. 

4.1.  Computing  nonlinear  internal  variables 

Reaction  flux:  Outputs  jDRA(z,t)  from  the  DRA  are  a  linearized 
approximation  to  the  true  j(z,t).  There  is  no  additional  correction  to 
this  variable. 

Overpotential:  We  find  that  outputs  t;dra (z,t)  from  the  DRA 
accurately  predict  the  variability  of  the  true  r\(z,t),  but  do  not  always 
predict  the  mean  of  ?j(z,t)  well.  The  average  value  for  p  across  the 
electrode,  based  on  the  linearized  model  of  Section  2,  is 


Potential  in  solid:  The  debiased  solid  potential  variable  0s(z,  t) 
is  defined  as  <j>s(z,  t)  —  (f>s(z.  t )  -  0S( 0,  t).  In  the  negative  electrode, 
4>s(0,t)  —  0  and  in  the  positive  electrode,  cj>s(0,t)  =  v(t),  the  cell 
voltage.  Therefore,  <j>s(z,t )  =  (ps(z,  t)  in  the  negative  electrode  and 
<j>s  (z,  t)  =  <j>s(z,  t)  +  v(t)  in  the  positive  electrode,  where  the  calcu¬ 
lation  of  v(t)  is  discussed  in  Section  4.2. 

Concentration  in  electrolyte:  The  electrolyte  concentration  is 
approximated  by  truncating  the  infinite-series  expansion  shown  in 
Eq.  (49).  The  equilibrium  value  is  added  to  the  appropriate  output  of 
the  state  space  model  to  arrive  at  the  true  value  of 
ce(x,t)  =  re(x,  t)  +  ce0. 

Solid-Electrolyte  potential  difference:  The  debiased 
solid— electrolyte  potential  difference  0s_e  is  defined  as 
<fis-e(z,t)  =  <l>s-e(z,t)  -  U0Cp(cSfi).  The  corresponding  transfer 
function  <£j_e(z,s)//app(s)  has  a  pole  at  s  =  0,  which  is  removed  to 
give  the  4>s_e(z.  s)//app(s)  transfer  function,  which  is  used  within 
the  DRA  to  produce  <j>s_e(z,t).  The  integrator  response  could  be 
added  back  manually,  as  in  0s_e(t)  =  0s-e(f)  +  (0s-e  )*i(t),  but 
better  performance  is  obtained  by  looking  deeper  at  what  is  actu¬ 
ally  happening. 

Recall  that  cSiavg(t)  =  (Sjee°)Xi(t)  +  cs0  and  note  that 

^resO  _  8U0cp|  x  gresO  J]lereforei  we  can  write 

0Cs-e  lcs,0 

0s -e(z,t)  =  0s-e(ziO+  ^Ef0cp(cSio)  +  |  qc°CP |  j  (Cs,avg  —  Cs,o)^  ’ 


linear (s)  =  ™ct  J  g^|/app(s)dz 
0 

^linear  (0  =  Q^fnf\  *aPpM: 


where  the  second  term  on  the  right-hand  side  is  equal  to  the  first 
two  terms  of  the  Taylor  series  expansion  of  liocp(cs,aVg)  around  the 
starting  average  concentration  cs,o.  Therefore,  instead  of  imple¬ 
menting  0s_e(z,s)//app(s),<we  find  that  we  achieve  more  accurate 
results  if  we  implement  0s_e(z,s)/I3pp(s )  and  then  compute 
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Me(M)  =  4>l-e(Z,t)  +  UocP(cs,mg(t)). 


Potential  in  electrolyte:  As  described  in  Section.  2.4,  the  elec¬ 
trolyte  potential  <|>e(x,t)  is  comprised  of  two  terms.  The  first  term 
[(|>e(x,t)]i  is  found  using  [0e(x,  t)],  which  is  an  output  from  the  DRA: 

[0e(x.t)]l  =  pe(*.0],+tfe(0,t). 

We  use  (t>e(0,t)  =  c[>s(0,t)— c(>s_e(0,t)  and  the  fact  that  c(>s(0,t)  —  0  at 
the  negative  current  collector  to  get 


0e(*.t)l  1  =  e(O-t)- 

The  second  term  in  the  c|>e(x,t)  equation,  [<pe(x,  t)]2,  is  given  by  Eq. 
(31).  The  values  for  ce(x,t)  are  outputs  of  the  state  space  model  with 
the  addition  of  the  equilibrium  value,  as  described  above.  The 
overall  electrolyte  potential  is  given  by 


MM)  =  [<Mx,t)]t  + 


2RT(1  -t°)  /Ce(x,m 

F  m\ce(0A)J- 


4.2.  Cell  voltage 

The  cell  voltage  is  the  difference  in  solid  potential  between  the 
two  current  collectors,  v(t)  =  <|>s(f.,t)-c|>s(0,t).  Using  the  over¬ 
potential  equation,  rj  =  c|>s-c|)e-Uocp-ERfiiin/  and 
<j>e(x,  t)  =  0e(x,  t)  -  0e(O,  t)  we  write 

v(t)  =  0s(L,t)-0s(O,t) 

=  7](L,  t)  +  ML,  t)  +  U0Cp(L,  t)  +  FRmm  pj(L,  t)  -  rj(0,  t) 


=  (t 7(L,  t)  -  77(0,  t))  +  ML,  t)  +  (U0cp(L,  t)  -  Uocp(0,  t)) 

+  t)  —  Ffllm.rjXO,  t))  . 

If  we  split  the  0e(L,  t)  term  into  its  two  components,  we  can 
write  the  output  as 

v(t )  =  F  (Rfllm,p j(L,  t)  -  Rfiim,rJ(0,  f))  +  [ML,  t)]  1  +  (y(L,  t) 

-i?(0,t))+  [0e(L,t)]2  +  (UOCp(cSle(I,t))  -l/ocp(cSle(0,t))). 

(57) 


At  first  glance,  it  appears  that  we  require  implementing 
a  minimum  of  nine  transfer  functions  to  be  able  to  compute  cell 
voltage:  j(L,t),  j(0,t),  [MU)],.  7,(0, t),  ce(L,t),  ce(0,t),  cs,e(L,t), 

and  cSie(0,t).  However,  the  linear  transfer  functions  can  be  grouped 
together.  Define 


VWs)  _  /  7(M_n  J(0,s)\  ,  [^L,s) j1 

/app(s)  V  fllm’P4pp(S)  mm’n/app(S);  +  fapp(s) 

/  i?(0,  s)\ 

V^app(s)  /app(S)/ 

If  this  transfer  function  is  implemented  in  the  DRA  to  produce 
a  single  output  viin(t),  then  we  can  write 


v(t)  =  Vlin(t)-%near,p(t)  +  t7nonlinear,p(i,f)  +  Wunear,n(t) 
-  Nonlinear, n(0,  t)  +  [^(L,  t)]  *  +  (U0cp  (cs,e(i,  t)) 
-l/ocp(Cs,e(0,t))). 


Table  1 

Cell  parameters  for  simulation. 


Symbol 

Units 

Negative 

electrode 

Separator 

Positive 

electrode 

L 

jim 

128 

76 

190 

R 

jim 

12.5 

- 

8.5 

A 

m2 

1 

1 

1 

a 

S  m 

1  100 

— 

3.8 

£S 

m3  m 

r3  0.471 

- 

0.297 

r3  0.357 

0.724 

0.444 

brug 

1.5 

- 

1.5 

c™ax 

n-3  26,390 

22,860 

ce,o 

n-3  2000 

2000 

2000 

0i,min 

- 

0.05 

- 

0.78 

0i.max 

- 

0.53 

- 

0.17 

Ds 

m2  s~ 

1  3.9  x  10~14 

— 

1.0  x  10-13 

De 

1  7.5  x  10"11 

7.5  x  10-11 

7.5  x  10-11 

t° 

k 

0.363 

0.363 

0.363 

1/2  m5/2  s-1  1.94  x  10-11 

- 

2.16  x  10-11 

0.5 

0.5 

Rfilm 

Qm2 

0.0 

- 

- 

We  compute  aeff  = 

:  ff£s,  =  K£Tas,  Df  =  De4rUS 

In  the 

electrolyte,  conductivity  is  a 

function  of 

concentration: 

«(Ce)  =  4. 

1253  x  1 
1.5094  x 

0-2  + 5.007  xl0-4ce- 4.7212 
10-'°c2  -  1.6018  x  10  14rf. 

x  lO-’c2 

For  the 

negative  electrode,  the  open-cir 

cuit  potential 

function  is: 

Uocp(0)  = 

-0.16  + 

1.32exp(—3.O0)+lO.Oexp(— 2000.00). 

For  the  positive  electrode,  the  open-circuit  potential  function  is:  t/ocp(®)  = 

4.19,829  +  0.0565661tanh(— 14.55460+8.60942)— 0.0275479  - 1 

[(0.998432-0)°  4924656 

— 1.901  llj  — 0.1 57123exp(—O.O473806)+O.81O239exp[— 40(0— 0.133875)]. 


This  output  equation  requires  implementing  a  minimum 
number  of  five  transfer  functions  to  be  able  to  compute  cell  voltage: 
Vlin(t),  Ce(I,t),  ce(0,t),  Cs,e(L,t),  and  CSre(0,t). 

5.  Simulation  and  results 

In  this  section,  we  use  simulation  to  demonstrate  the  perfor¬ 
mance  of  the  reduced-order  model  as  compared  to  the  full-order 
PDEs  describing  the  porous-electrode  model.  The  parameters  of 
the  cell  being  simulated  are  those  published  by  Doyle  et  al.  [15],  and 
are  listed  in  Table  1.  The  cell  input  current  for  the  simulation  is 
based  on  the  Environmental  Protection  Agency’s  (EPA)  Urban 
Dynamometer  Driving  Schedule  (UDDS),  which  was  developed  by 
the  EPA  to  represent  city  driving  conditions  for  light  duty  vehicles. 
Fig.  2  shows  the  normalized  UDDS  profile  of  current  versus  time,  as 
computed  for  a  medium-size  electric  vehicle.  For  our  simulations 
we  multiplied  this  profile  by  two,  to  achieve  a  maximum  absolute 
rate  of  2C,  which  for  this  cell  corresponds  to  41  A. 


(58) 


Fig.  2.  Normalized  profile  of  cell  current  ’ 


for  the  UDDS  cycle. 
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We  compare  our  ROM  results  to  simulations  of  the  pseudo-2D 
porous-electrode  model,  which  we  implemented  in  COMSOL 
Multiphysics  version  4.2a.  We  implemented  two  different  versions 
of  the  porous-electrode  model:  (1)  a  “rigorous  PDE”  version,  which 
implements  Eqs.  (1)— (6)  as  listed,  and  (2)  a  “linearized  PDE” 
version,  which  linearizes  the  PDEs  of  Eqs.  (l)-(6)  to  exactly  match 
Assumptions  1  and  2.  Specifically,  the  linearized  PDE  model  is 
created  from  the  rigorous  PDE  model  with  the  following  modifi¬ 
cations:  (1)  the  Deep  functions  are  linearized  about  the  initial  cell 


Fig.  4.  Eigenfunctions  used  in  the  computation  of  electrolyte  concentration. 


SOC,  (2)  the  second  term  of  Eq.  (15)  is  dropped,  (3)  the  value  of  k 
remains  constant  throughout  the  simulation,  at  the  value  corre¬ 
sponding  to  the  initial  electrolyte  concentration,  and  (4)  we  use  the 
linearized  Butler— Volmer  equation  with  the  jo  term  constant 
throughout  the  simulation.  Because  these  same  assumptions  are 
made  to  generate  the  ROM,  we  expect  that  ROM  solution  will 
converge  to  the  linearized  PDE  model  as  we  increase  the  order  of 
the  ROM. 

In  all  cases,  the  simulation  is  initialized  with  a  cell  SOC  of  60%. 
We  compute  DRA  outputs  and  PDE  outputs  for  all  seven  cell  vari¬ 
ables  (0s-e,  y,j,  cs,  e,  </> s,  0e,  and  ce)  at  four  different  spatial  locations 
across  the  cell  (at  both  current  collectors,  and  at  both  electrode/ 
separator  interfaces).  The  reduced-order  state-space  model  there¬ 
fore  has  one  input  (cell  current)  and  28  outputs. 

Fig.  3  shows  the  frequency  responses  for  two  of  the  continuous¬ 
time  transfer  functions.  The  reaction  flux  is  shown  in  Frame  (a) 
which  is  a  plot  of  Eq.  (21)  at  z  locations  ranging  from  0  to  1  at  0.1 
increments.  The  surface  concentration,  with  the  pole  at  s  =  0 
removed  (Eq.  (51)),  is  plotted  in  Frame  (b)  at  the  same  z  locations. 
We  see  that  although  the  transfer  functions  are  mathematically 
complex,  the  frequency  response  of  the  variables  is  actually  quite 
simple.  This  is  why  a  reduced-order  model  can  accurately  approx¬ 
imate  the  rigorous  PDE  model. 

When  approximating  the  electrolyte  concentration,  we  use 
eigenvalue/eigenfunction  pairs  4  and  W(X'M  for  n  =  0,...,5.  The 
“zeroth”  eigenvalue/eigenfunction  simply  gives  the  equilibrium 


Ordered  singular  values  Cell  voltage  RMS  error  for  block  Hankel  size  =  8000 


o  ROM  vereu"  nneari2ed>PDE| 

I6 

o  5 

RMS  eri 

‘  .  ■ . 

123456789  10 

Singular  value 


23456789  10  11 

System  order 


(a) 


(b) 


Fig.  5.  Singular  values  in  DRA  step  3  for  a  Hankel  matrix  size  of  8000,  and  RMS  voltage  prediction 


:  for  different  system 
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Fig.  6.  RMS  voltage  prediction  errors  for  different  Hankel  matrix  sizes  in  DRA  step  3. 


concentration  (ce,o  =  2000  mol  m  3)  and  is  therefore  not  included 
in  the  state-space  model  explicitly.  Fig.  4  plots  eigenfunctions 
W(x\Xk)  for  n  =  1....5. 

Four  tuning  parameters  must  be  set  to  generate  a  cell  model 
using  the  DRA:  the  high-rate  sampling  frequency  Fs  for  Step  1  of  the 
DRA,  the  size  of  the  Hankel  matrix  for  Step  3  of  the  DRA,  the 
number  of  states  in  the  reduced-order  system,  and  the  sampling 
period  Ts  of  the  computed  discrete-time  reduced-order  model.  For 
the  simulation  runs  presented  in  this  paper,  the  high-rate  sampling 
frequency  was  chosen  as  Fs  =  128  Hz  and  the  sampling  period  of  the 
computed  discrete-time  reduced-order  model  was  chosen  as 
Ts  =  Is.  The  computational  complexity  required  to  run  the  DRA  and 
the  accuracy  of  the  resulting  reduced-order  model  predictions  are 
not  sensitive  to  either  of  these  parameters. 

We  use  the  ordered  singular  values  of  the  system  Hankel  matrix 
in  Step  3  of  the  DRA  to  indicate  the  relative  importance  of  each 
state,  to  assist  in  determining  the  best  tradeoff  between  final  ROM 
complexity  and  fidelity.  Fig.  5(a)  shows  the  first  ten  singular  values. 
The  tradeoff  between  accuracy  and  complexity  is  application 
dependent  but  it  is  evident  that  there  is  no  need  to  include  any 
more  than  the  sixth  singular  value  (which  is  about  1%  of  the  first 
singular  value).  The  root-mean-squared  (RMS)  error  between  the 
PDE  simulation  and  the  ROM  simulation  of  cell  voltage  as  a  function 
of  the  system  order  is  shown  in  Fig.  5(b).4  We  see  that  the  RMS  error 


4  Note  the  system  order  in  this  plot  includes  the  addition  of  the  augmented 
integrator  term. 


between  the  ROM  and  the  linearized  PDE  model  approaches  zero, 
and  that  the  RMS  error  between  the  ROM  and  the  rigorous  PDE 
model  approaches  a  constant  value.  In  both  cases,  we  achieve  little 
improvement  in  RMS  error  after  about  the  fifth  state.  Thus,  from 
now  on,  we  use  a  ROM  having  exactly  five  state  variables  in  its  state 
vector  x(t). 

The  final  DRA  tuning  parameter  is  the  length  of  the  discrete¬ 
time  pulse  response  included  in  the  Hankel  matrix.  This  value 
must  be  large  enough  to  capture  the  extent  of  the  longest  pulse 
response  (slowest  system  dynamics)  which  turn  out  to  be  that  of 
the  surface  concentration  of  lithium  in  the  positive  electrode.  Fig.  6 
shows  that  the  cell  voltage  RMS  error  is  fairly  invariant  to  the  block 
Hankel  matrix  size,  although  we  find  that  the  fidelity  of  individual 
transfer  functions  can  be  sensitive  to  this  value,  and  larger  lengths 
usually  help. 

Based  on  these  results,  we  generate  a  ROM  having  five  states 
using  a  Hankel  matrix  block  size  of  8000.  Fig.  7  compares  the 
reaction  flux  and  the  solid  surface  concentration  of  the  ROM  to  the 
linearized  PDE  model.  Although  the  figure  shows  results  for  only 
two  of  the  28  outputs,  the  other  26  match  similarly  well.  This 
indicates  that  we  can  reduce  the  infinite-order  system  to  a  fifth- 
order  system  with  high  accuracy. 

Results  of  comparing  this  ROM  to  the  rigorous  PDE  model  at 
the  electrode/separator  interfaces  are  presented  in  Fig.  8.  The 
electrochemical  reactions  at  the  electrode/separator  interface 
are  more  dynamic  than  at  the  current  collectors  and  are  there¬ 
fore  more  difficult  for  the  ROM  to  match.  Of  particular  interest 
are  the  results  for  the  reaction  flux  (Frames  (a)  and  (b))  where 
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Fig.  8.  Comparison  of  the  fifth-c 


ROM 
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c  at  negative  current  collector  c  at  positive  current  collector 


the  ROM  does  not  match  the  rigorous  PDE  model  especially  well 
at  highly  dynamic  times  (e.g.,  near  t  =  170  s).  We  attribute  this 
error  to  the  assumptions  made  in  linearizing  the  full-order 
model.  The  rigorous  PDE  model  predicts  a  higher  dynamic 
response  than  does  the  linearized  PDE  model,  which  can  be  seen 
by  comparing  Fig.  8(b)  and  Fig.  7(a).  The  same  effect  is  seen  with 
the  surface  lithium  concentration  at  the  positive  electrode/ 
separator  boundary  which  is  shown  in  Fig.  8(h).  Again  the  ROM 
predicts  the  linearized  PDE  response  extremely  well  as  shown  in 
Fig.  7(b),  but  the  linearized  PDE  model  and  the  rigorous  PDE 
model  differ. 

Increasing  the  model  order  or  increasing  the  size  of  the  Hankel 
matrix  does  not  improve  modeling  of  j  or  cs,e  because  the  ROM 
already  matches  the  linearized  PDE  model.  However,  unlike  j  or  cs,e, 


the  results  for  the  electrolyte  concentration  can  be  improved  by 
increasing  the  system  order.  The  infinite-series  expansion  used  to 
represent  the  electrolyte  concentration  is  truncated  in  the  ROM. 
Increasing  the  number  of  states  results  in  slightly  better  estimation 
of  ce.  Fig.  9(a— d)  shows  the  comparison  between  the  fifth-order 
ROM  and  the  rigorous  PDE  model.  Fig.  9(e,f)  shows  the  concen¬ 
tration  at  the  separator  interface  with  an  eighth-order  ROM. 

The  output  cell  voltage  for  the  ROM  is  calculated  from 
the  electrochemical  variables  using  Eq.  (57).  The  results  of  the  fifth- 
order  ROM  and  the  rigorous  PDE  model  are  shown  in  Fig.  10.  The 
cell  voltage  RMS  error  between  the  rigorous  PDE  model  and  the 
ROM  is  1.14  mV.  Frame  (b)  shows  a  detailed  view  of  the  simulation 
during  the  period  of  the  greatest  mismatch.  Frame  (c)  illustrates  the 
difference  between  the  corrected  and  uncorrected  77  value,  as 
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Eq.  (58)  can  be  used  to  turn  these  five  linear  outputs  into 
Table2  a  nonlinear  voltage  estimate.  Simulations  using  this  minimum- 

rom  for  five  states,  to  compute  output  voltage  only.  output  fifth-order  reduced-order  model  match  those  shown  in  Fig.  10. 


0.7568  0  0  0 

0  0.9840  0  0 

A  =  0  0  0.9969  0 

0  0  0  0.9998 

0  0  0  0 


—6.459 x  10-5  — 7.214x  10“6  -1.280xl0-6  -3.572xl0-? 
3 .489 x  10-2  2.861  x  10-2  -8.019xl0-2  -7.194xl0-4 

C=  -1.287x10-'  3.944x  10“2  1.149x10-'  -1.171xl0-3 

5.569x10“'  3.514xl0-2  9.591  xl0“3  -9.381  xl0“2 

-1.300x10+°  -2.369x10“'  -3.840xlQ-2  1.692xlQ-3 


discussed  in  Section.  4,  over  the  same  period  and  Frame  (d)  shows 
this  impact  on  the  cell  voltage. 

The  fifth-order  ROM  used  in  these  simulations  has  28  linear 
outputs.  In  general,  the  cell  voltage  calculation  itself  requires  only 
five  linear  outputs.  Table  2  shows  the  A,  B,  C,  and  D  matrices  for 
a  minimum-output  fifth-order  reduced-order  model  that  can  be 
used  to  calculate  cell  output  voltage.  The  input  is  iapp(t).  and  the 
output  is 


"lin(t) 
Ce(L,  t) 
Ce(0,t) 
Cs,e(L,t) 
Cs  e1  0.1) 


6.  Conclusions 

In  this  paper  we  have  developed  a  one-dimensional  physics- 
based  reduced-order  model  for  lithium-ion  cells.  To  do  so,  we  first 
augmented  the  set  of  electrochemical-variable  transfer  functions 
developed  by  Smith  et  al.  [8,9]  with  transfer  functions  for  the  solid 
and  electrolyte  potentials  c>s  and  4>e.  and  the  electrolyte  concen¬ 
tration  ce.  The  discrete-time  realization  algorithm  is  used  to 
produce  an  optimal  reduced-order  discrete-time  state-space  model 
from  these  transfer  functions,  overcoming  the  limitations  of 
nonlinear  optimization  approaches.  Nonlinear  corrections  to  the 
linearized  model  are  also  given  to  improve  accuracy  of  the  method. 

The  final  state-space  ROM  used  in  Section.  5  requires  134  msec  to 
simulate  the  1500  s  UDDS  profile  on  a  desktop  computer,  which  is 
about  0.09  msec  sample-1.  In  contrast,  the  rigorous  model,  which  is 
run  using  COMSOL,  requires  almost  13  min  (or  510  msec  sample-1) 
for  the  same  simulation  on  the  same  desktop  computer.  This  is 
a  speedup  of  over  5000:1.  Thus,  we  conclude  that  the  simplicity  of 
the  ROM  makes  it  feasible  for  use  in  real-time  control  applications 
such  as  in  electric-vehicle  battery-management  systems. 

Note  that  a  change  in  parameter  values,  such  as  happens  during 
cell  aging,  does  not  change  the  transfer-function  equations  of  the 
model  developed  in  Section.  2;  however,  the  state-space  reduced- 
order  model  will  need  to  be  regenerated  from  these  transfer 
functions  using  the  DRA  discussed  in  Section.  3.  The  DRA  requires 
about  11  min  on  the  same  desktop  computer  to  generate  the 
reduced-order  model  from  the  transfer  functions.  (Again,  this 
process  must  be  performed  only  ( 1 )  once  to  create  the  original  ROM 
from  the  original  set  of  parameters,  and  (2)  whenever  the 
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parameter  values  are  deemed  to  have  changed  significantly  from 
their  original  values — perhaps  once  per  month  of  regular  operation, 
and  when  it  can  be  performed  during  an  interval  when  the  BMS  is 
otherwise  essentially  idle,  as  might  be  the  case  when  the  battery 
pack  is  being  charged.)  We  note  that  the  extended  time  required  to 
calculate  the  ROM  when  parameter  values  are  modified  makes  it 
impractical  for  use  in  system-identification  applications. 

Finally,  the  results  show  that  a  fifth-order  ROM  is  able  to  track  the 
highly  dynamic  UDDS  profile  of  current  versus  time  with  an  output 
cell  voltage  root-mean-squared  error  of  less  than  1.2  mV  for  the 
specific  cell  parameters  considered.  This  ROM  also  tracks  all  of  the 
cell’s  electrochemical  states  at  any  position  across  its  cross  section 
with  a  high  degree  of  accuracy.  The  greatest  source  of  mismatch 
between  the  ROM  and  PDE  simulation  predictions  arises  as  a  result 
of  the  linearization  assumptions  made  when  creating  the  transfer 
functions  of  the  model.  As  the  cell  state  moves  away  from  from  the 
initial  linearization  set-point,  prediction  errors  also  increase.  In 
a  future  paper  we  will  discuss  a  method  to  overcome  this  limitation. 
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Nomenclature5 

A:  surface  area  of  the  electrode,  m2 

A:  state  transition  matrix  of  the  state-space  model 


as:  specific  surface  area  of  the  porous  electrode,  m2  m-3 
B :  input  matrix  of  the  state-space  model 
C:  output  matrix  of  the  state-space  model 

c:  concentration  of  lithium  in  phase  indicated  by  subscript,  mol  m-3 
ce,o.’  steady-state  concentration  of  lithium  in  the  electrolyte  phase,  mol  m 
c*  k:  kth  generalized  Fourier  coefficient  of  ce(x,t) 

n  lithium  concentration  in  an  electrode  particle, ^mol  m-3 

cs,ovg.'  average  concentration  of  lithium  in  an  electrode  particle,  mol  m-3 

cs  e:  surface  concentration  of  lithium  in  a  spherical  electrode  particle,  mol 

Cs,e.‘  debiased  version  of  cSte,  cs,e  =  cs,e  -  cs  0.  mol  m-3 

D:  input— output  coupling  matrix  of  the  state-space  model 

of :  effective  electrolyte  diffusivity,  m2  s_1 

Ds:  solid  diffusivity,  m2  s_1 

f±:  mean  molar  activity  coefficient 

F:  Faraday's  constant,  96  487  C  mol-1 

Fs:  sampling  frequency  of  the  transfer  functions  in  DRA  Step  1 
1:  applied  current  density,  A  m-2 
iap P:  applied  cell  current,  iapp  =  IA,  A 

j:  reaction  flux,  mol  m-2  s_1 

fk:  kth  generalized  Fourier  coefficient  of  J(x,t) 

k:  rate  constant  for  the  electrochemical  reaction,  mol“_1m4_3“s_1 

L:  (without  subscript)  length  of  the  cell,  m 

L:  (with  subscript)  length  of  region  of  cell,  m 

rs:  sampling  period  of  final  discrete-time  state  space  system,  s 

r:  radial  coordinate,  m 

reso:  residues  of  poles  of  transfer  functions  at  s  =  0 
R:  universal  gas  constant,  8.31451  J  mol-1  K_1 
Ra:  charge  transfer  resistance,  Q  m2 
Rfilm:  film  resistance,  Q  m2 
Rs:  particle  radius,  m 

Rm:  sum  of  the  charge  transfer  resistance  and  the  film  resistance,  £1  m2 

T:  temperature,  K 

t°:  transference  number 


UoqJ:  open  circuit  potential,  V 
v:  cell  voltage,  V 

x:  ID  linear  coordinate  across  the  cell,  m 
x:  state  of  state-space  model 
y:  linear  output  of  state-space  model 
z:  unitless  dimension  across  electrode 


a:  charge-transfer  coefficient 

fi:  parameter  of  Jacobsen-West  transfer  function:  (J  =  Rss/Ds 
r:  volume  fraction  of  phase  indicated  by  subscript 
<ji:  potential  of  the  phase  indicated  by  subscript,  V 
<j>:  debiased  potential  of  the  phase  indicated  by  subscript,  V 
0s_e:  potential  difference  between  solid  and  electrolyte  phases,  V 
l<j>eh-  linear  part  of  c[>e,  V 
[(Jeh-  nonlinear  part  of  (|)e,  V 
rf.  local  overpotential,  V 
ksS:  effective  electrolyte  conductivity,  S  m_1 
fa:  fcth  eigenvalue  of  the  electrolyte  concentration  function 
v:  non-spatial-dependent  ODE  coefficient,  dimensionless 
V:  eigenfunction  of  the  electrolyte  concentration  function 
effective  solid  conductivity,  S  m_1 

Subscript/superscript 

e:  pertaining  to  the  electrolyte  phase 

m:  pertaining  to  the  separator  ("middle”  of  cell) 

k:  pertaining  to  the  negative  electrode 

p:  pertaining  to  the  positive  electrode 

s:  pertaining  to  the  solid  phase 


5  Unless  otherwise  specified,  lowercase  symbols  correspond  to  constants  or 
time-domain  quantities;  uppercase  symbols  correspond  to  Laplace  frequency- 


